/**********************************************************************

  Audacity: A Digital Audio Editor

  FreqWindow.cpp

  Dominic Mazzoni

  *******************************************************************//**

  \class FreqWindow
  \brief Displays a spectrum plot of the waveform.  Has options for
  selecting parameters of the plot.

  Has a feature that finds peaks and reports their value as you move
  the mouse around.

  *//****************************************************************//**

  \class FreqPlot
  \brief Works with FreqWindow to dsplay a spectrum plot of the waveform.
  This class actually does the graph display.

  Has a feature that finds peaks and reports their value as you move
  the mouse around.

  *//*******************************************************************/

/*
  Salvo Ventura - November 2006
  Extended range check for additional FFT windows
  */

#ifdef __BORLANDC__
#pragma hdrstop
#endif

#include <math.h>
#include <string.h>
#include <memory.h>
#include "FFT.h"
#include "qfreqalg.h"

#ifndef min
#define min(a,b)            (((a) < (b)) ? (a) : (b))
#endif

// FreqWindow
namespace AudioUtils {

	QFreqAlg::QFreqAlg()
		: mData(NULL)
		, mProcessed(NULL)
		, mWindowSize(1024)
	{
		mRate = 0;
		mDataLen = 0;

		mAlg = 0;
		mFunc = 0;

		dBRange = 120;
		mLogAxis = false;
	}

	QFreqAlg::~QFreqAlg()
	{
		if (mData)
			delete[] mData;

		if (mProcessed)
			delete[] mProcessed;

		DeinitFFT();
	}

	void QFreqAlg::setRate(double rate)
	{
		mRate = rate;
	}

	void QFreqAlg::setAlg(int alg)
	{
		mAlg = alg;
	}

	void QFreqAlg::setFunc(int func)
	{
		mFunc = func;
	}

	void QFreqAlg::setWindowSize(int size)
	{
		mWindowSize = size;
	}

	void QFreqAlg::setAxis(int axis)
	{
		mLogAxis = axis ? true : false;
	}

	void QFreqAlg::setData(float* data, int len)
	{
		if (mData) {
			delete[] mData;
		}

		mDataLen = len;
		mData = new float[mDataLen];

		memset(mData, 0, mDataLen * sizeof(float));
		memcpy(mData, data, mDataLen * sizeof(float));
	}

	void QFreqAlg::setData(short* data, int len)
	{
		double fscale = 1.0 / 32768.0;
		float* fdata = new float[len];
		for (int i = 0; i < len; i++)
		{
			fdata[i] = (float)((double)data[i] * fscale);
		}

		setData(fdata, len);
		delete[] fdata;
	}

	float QFreqAlg::cubicInterpolate(float y0, float y1, float y2, float y3, float x)
	{
		float a, b, c, d;

		a = y0 / -6.0 + y1 / 2.0 - y2 / 2.0 + y3 / 6.0;
		b = y0 - 5.0 * y1 / 2.0 + 2.0 * y2 - y3 / 2.0;
		c = -11.0 * y0 / 6.0 + 3.0 * y1 - 3.0 * y2 / 2.0 + y3 / 3.0;
		d = y0;

		float xx = x * x;
		float xxx = xx * x;

		return (a * xxx + b * xx + c * x + d);
	}

	float QFreqAlg::cubicMaximize(float y0, float y1, float y2, float y3, float * max)
	{
		// Find coefficients of cubic

		float a, b, c, d;

		a = y0 / -6.0 + y1 / 2.0 - y2 / 2.0 + y3 / 6.0;
		b = y0 - 5.0 * y1 / 2.0 + 2.0 * y2 - y3 / 2.0;
		c = -11.0 * y0 / 6.0 + 3.0 * y1 - 3.0 * y2 / 2.0 + y3 / 3.0;
		d = y0;

		// Take derivative

		float da, db, dc;

		da = 3 * a;
		db = 2 * b;
		dc = c;

		// Find zeroes of derivative using quadratic equation

		float discriminant = db * db - 4 * da * dc;
		if (discriminant < 0.0)
			return float(-1.0);              // error

		float x1 = (-db + sqrt(discriminant)) / (2 * da);
		float x2 = (-db - sqrt(discriminant)) / (2 * da);

		// The one which corresponds to a local _maximum_ in the
		// cubic is the one we want - the one with a negative
		// second derivative

		float dda = 2 * da;
		float ddb = db;

		if (dda * x1 + ddb < 0)
		{
			*max = a*x1*x1*x1 + b*x1*x1 + c*x1 + d;
			return x1;
		}
		else
		{
			*max = a*x2*x2*x2 + b*x2*x2 + c*x2 + d;
			return x2;
		}
	}

	float QFreqAlg::getProcessedValue(float freq0, float freq1)
	{
		int alg = mAlg;

		float bin0, bin1, binwidth;

		if (alg == 0) {
			bin0 = freq0 * mWindowSize / mRate;
			bin1 = freq1 * mWindowSize / mRate;
		}
		else {
			bin0 = freq0 * mRate;
			bin1 = freq1 * mRate;
		}
		binwidth = bin1 - bin0;

		float value = float(0.0);

		if (binwidth < 1.0) {
			float binmid = (bin0 + bin1) / 2.0;
			int ibin = int(binmid) - 1;
			if (ibin < 1)
				ibin = 1;
			if (ibin >= mProcessedSize - 3)
				ibin = mProcessedSize - 4;

			value = cubicInterpolate(mProcessed[ibin],
				mProcessed[ibin + 1],
				mProcessed[ibin + 2],
				mProcessed[ibin + 3], binmid - ibin);

		}
		else {
			if (int(bin1) > int(bin0))
				value += mProcessed[int(bin0)] * (int(bin0) + 1 - bin0);
			bin0 = 1 + int(bin0);
			while (bin0 < int(bin1)) {
				value += mProcessed[int(bin0)];
				bin0 += 1.0;
			}
			value += mProcessed[int(bin1)] * (bin1 - int(bin1));

			value /= binwidth;
		}

		return value;
	}

	void QFreqAlg::recalc()
	{
		if (mProcessed)
			delete[] mProcessed;
		mProcessed = NULL;

		if (!mData) {
			return;
		}

		int f = 10;
		int alg = mAlg;
		int windowFunc = mFunc;
		long windowSize = mWindowSize;

		if (!(windowSize >= 32 && windowSize <= 65536 &&
			alg >= 0 && alg <= 4 && windowFunc >= 0 && windowFunc < f)) {
			return;
		}

		if (mDataLen < mWindowSize) {
			return;
		}

		mProcessed = new float[mWindowSize];

		int i;
		for (i = 0; i < mWindowSize; i++)
			mProcessed[i] = float(0.0);
		int half = mWindowSize / 2;

		float *in = new float[mWindowSize];
		float *in2 = new float[mWindowSize];
		float *out = new float[mWindowSize];
		float *out2 = new float[mWindowSize];
		float *win = new float[mWindowSize];

		// initialize the window
		for (int i = 0; i < mWindowSize; i++)
			win[i] = 1.0;
		WindowFunc(windowFunc, mWindowSize, win);
		// Scale window such that an amplitude of 1.0 in the time domain
		// shows an amplitude of 0dB in the frequency domain
		double wss = 0;
		for (int i = 0; i < mWindowSize; i++)
			wss += win[i];
		if (wss > 0)
			wss = 4.0 / (wss*wss);
		else
			wss = 1.0;

		int start = 0;
		int windows = 0;
		while (start + mWindowSize <= mDataLen) {
			for (i = 0; i < mWindowSize; i++)
				in[i] = win[i] * mData[start + i];

			switch (alg) {
			case 0:                  // Spectrum
				PowerSpectrum(mWindowSize, in, out);

				for (i = 0; i < half; i++)
					mProcessed[i] += out[i];
				break;

			case 1:
			case 2:
			case 3:   // Autocorrelation, Cuberoot AC or Enhanced AC

				// Take FFT
#ifdef EXPERIMENTAL_USE_REALFFTF
				RealFFT(mWindowSize, in, out, out2);
#else
				FFT(mWindowSize, false, in, NULL, out, out2);
#endif
				// Compute power
				for (i = 0; i < mWindowSize; i++)
					in[i] = (out[i] * out[i]) + (out2[i] * out2[i]);

				if (alg == 1) {
					for (i = 0; i < mWindowSize; i++)
						in[i] = sqrt(in[i]);
				}
				if (alg == 2 || alg == 3) {
					// Tolonen and Karjalainen recommend taking the cube root
					// of the power, instead of the square root

					for (i = 0; i < mWindowSize; i++)
						in[i] = pow(in[i], 1.0f / 3.0f);
				}
				// Take FFT
#ifdef EXPERIMENTAL_USE_REALFFTF
				RealFFT(mWindowSize, in, out, out2);
#else
				FFT(mWindowSize, false, in, NULL, out, out2);
#endif

				// Take real part of result
				for (i = 0; i < half; i++)
					mProcessed[i] += out[i];
				break;

			case 4:                  // Cepstrum
#ifdef EXPERIMENTAL_USE_REALFFTF
				RealFFT(mWindowSize, in, out, out2);
#else
				FFT(mWindowSize, false, in, NULL, out, out2);
#endif

				// Compute log power
				// Set a sane lower limit assuming maximum time amplitude of 1.0
				float power;
				float minpower = 1e-20*mWindowSize*mWindowSize;
				for (i = 0; i < mWindowSize; i++)
				{
					power = (out[i] * out[i]) + (out2[i] * out2[i]);
					if (power < minpower)
						in[i] = log(minpower);
					else
						in[i] = log(power);
				}
				// Take IFFT
#ifdef EXPERIMENTAL_USE_REALFFTF
				InverseRealFFT(mWindowSize, in, NULL, out);
#else
				FFT(mWindowSize, true, in, NULL, out, out2);
#endif

				// Take real part of result
				for (i = 0; i < half; i++)
					mProcessed[i] += out[i];

				break;
			}                         //switch

			start += half;
			windows++;
		}

		switch (alg) {
			double scale;
		case 0:                     // Spectrum
			// Convert to decibels
			mYMin = 1000000.;
			mYMax = -1000000.;
			scale = wss / (double)windows;
			for (i = 0; i < half; i++)
			{
				mProcessed[i] = 10 * log10(mProcessed[i] * scale);
				if (mProcessed[i] > mYMax)
					mYMax = mProcessed[i];
				else if (mProcessed[i] < mYMin)
					mYMin = mProcessed[i];
			}
			if (mYMin < -dBRange)
				mYMin = -dBRange;
			if (mYMax <= -dBRange)
				mYMax = -dBRange + 10.; // it's all out of range, but show a scale.
			else
				mYMax += .5;

			mProcessedSize = half;
			mYStep = 10;
			break;

		case 1:                     // Standard Autocorrelation
		case 2:                     // Cuberoot Autocorrelation
			for (i = 0; i < half; i++)
				mProcessed[i] = mProcessed[i] / windows;

			// Find min/max
			mYMin = mProcessed[0];
			mYMax = mProcessed[0];
			for (i = 1; i < half; i++)
				if (mProcessed[i] > mYMax)
					mYMax = mProcessed[i];
				else if (mProcessed[i] < mYMin)
					mYMin = mProcessed[i];

			mYStep = 1;

			mProcessedSize = half;
			break;

		case 3:                     // Enhanced Autocorrelation
			for (i = 0; i < half; i++)
				mProcessed[i] = mProcessed[i] / windows;

			// Peak Pruning as described by Tolonen and Karjalainen, 2000

			// Clip at zero, copy to temp array
			for (i = 0; i < half; i++) {
				if (mProcessed[i] < 0.0)
					mProcessed[i] = float(0.0);
				out[i] = mProcessed[i];
			}

			// Subtract a time-doubled signal (linearly interp.) from the original
			// (clipped) signal
			for (i = 0; i < half; i++)
				if ((i % 2) == 0)
					mProcessed[i] -= out[i / 2];
				else
					mProcessed[i] -= ((out[i / 2] + out[i / 2 + 1]) / 2);

			// Clip at zero again
			for (i = 0; i < half; i++)
				if (mProcessed[i] < 0.0)
					mProcessed[i] = float(0.0);

			// Find new min/max
			mYMin = mProcessed[0];
			mYMax = mProcessed[0];
			for (i = 1; i < half; i++)
				if (mProcessed[i] > mYMax)
					mYMax = mProcessed[i];
				else if (mProcessed[i] < mYMin)
					mYMin = mProcessed[i];

			mYStep = 1;

			mProcessedSize = half;
			break;

		case 4:                     // Cepstrum
			for (i = 0; i < half; i++)
				mProcessed[i] = mProcessed[i] / windows;

			// Find min/max, ignoring first and last few values
			int ignore = 4;
			mYMin = mProcessed[ignore];
			mYMax = mProcessed[ignore];
			for (i = ignore + 1; i < half - ignore; i++)
				if (mProcessed[i] > mYMax)
					mYMax = mProcessed[i];
				else if (mProcessed[i] < mYMin)
					mYMin = mProcessed[i];

			mYStep = 1;

			mProcessedSize = half;
			break;
		}

		delete[]in;
		delete[]in2;
		delete[]out;
		delete[]out2;
		delete[]win;
	}

	int QFreqAlg::getFreqDbSize() const
	{
		return mProcessedSize;
	}

	const float* QFreqAlg::getFreqDb() const
	{
		return mProcessed;
	}

	void QFreqAlg::getFreqDb(float* buffer, int len)
	{
		int nCount = min(mProcessedSize, len);
		if (nCount > 0)
		{
			memcpy(buffer, mProcessed, nCount * sizeof(float));
		}
	}

	float QFreqAlg::getYMin() const
	{
		return mYMin;
	}

	float QFreqAlg::getYMax() const
	{
		return mYMax;
	}
}